set.seed(1)
#if (!requireNamespace("BiocManager", quietly = TRUE))   
#  install.packages("BiocManager",repos = "http://cran.us.r-project.org")  
#BiocManager::install("fgsea")
#BiocManager::install("limma")
#BiocManager::install("Biobase")
#install.packages("ggplot2")
#install.packages("igraph")
library(limma)
library(Biobase)
library(ggplot2)
library(igraph)
library(data.table)
library(fgsea)
library(readr)
library(biomaRt)
library(edgeR)

##############################
# Load KEGG pathways
pathways <- gmtPathways("/home/ubuntu/c2.cp.kegg.v7.1.symbols.gmt")
sig_path_limma = read.delim("/home/esaha/BONOBO/plots/GSE19783_miRNA_sigPathways_limma_lumAB.txt")[,1]
pathgenes = unique(unlist(pathways[which(names(pathways) %in% sig_path_limma)]))


#bonobo = data.frame(fread("/home/esaha/BONOBO/network/miRNA_mRNA_breastCancer/lumAB_edges_immueGenes.txt"))
bonobo = data.frame(fread("/home/esaha/BONOBO/network/miRNA_mRNA_breastCancer/GSE19783_miRNAmRNA_edges.txt"))

edges = bonobo$V1
genes = unlist(lapply(strsplit(edges, split = ":"), function(x){x[2]}))

indegree = data.frame(bonobo[which(genes %in% pathgenes),-1])
rownames(indegree) = edges[which(genes %in% pathgenes)]

# Get phenotypic data
phenotype = read.csv("~/BONOBO/data/miRNA_mRNA_breastCancer/GSE19783_phenotype.txt", sep="")

subtypes = phenotype$breast.cancer.subtype.ch1[match(colnames(indegree), rownames(phenotype))]
indegree = indegree[,which(subtypes %in% c("Lum A", "Lum B"))]
subtypes = subtypes[which(subtypes %in% c("Lum A", "Lum B"))]
subtypes = factor(subtypes)

# limma model: rank by intercept for each subtype for GSEA (don't include population inteecept)
#### Design Matrix

design = model.matrix(~ subtypes)
fit <- lmFit(indegree, design)
fit <- eBayes(fit)

tb = topTable(fit,coef="subtypesLum B",number=Inf)
head(tb)

sigEdges = rownames(tb)[which(tb$P.Val < 0.00005)]
weights = tb$t[which(tb$P.Val < 0.00005)]

sig_genes = unlist(lapply(strsplit(sigEdges, split = ":"), function(x){x[2]}))
sig_miRNA = unlist(lapply(strsplit(sigEdges, split = ":"), function(x){x[1]}))

unique(sig_genes)
unique(sig_miRNA)

# Keep miRNA with degree > 5 or gene with degree > 1
keep_miRNA = names(table(sig_miRNA))[as.vector(table(sig_miRNA)) > 4]
# keep_gene = names(table(sig_genes))[as.vector(table(sig_genes)) > 1]

sig_genes = sig_genes[which(sig_miRNA %in% keep_miRNA)]
sig_miRNA = sig_miRNA[which(sig_miRNA %in% keep_miRNA)]
weights = weights[which(sig_miRNA %in% keep_miRNA)]

diffNet = data.frame(weights)
colnames(diffNet) = "force"
diffNet$from = sig_miRNA
diffNet$to = sig_genes
diffNet$arrows = "to" 
head(diffNet)

# Plot graph
library(visNetwork)

# Network Plot
nodes <- data.frame(id = unique(as.vector(as.matrix(diffNet[,c(2,3)]))) , 
                    label=unique(as.vector(as.matrix(diffNet[,c(2,3)]))))
nodes$group <- ifelse(nodes$id %in% diffNet$from, "miRNA", "gene")

# Color male edges blue and female edges red
edge_col = rep("red", nrow(diffNet))
edge_col[which(sign(diffNet[,1]) == -1)] = "blue"

diffNet$color = edge_col

net <- visNetwork(nodes, diffNet, width = "100%", 
                  main = "Most Differential Edges (Lum B - Lum A)") 
net <- visGroups(net, groupname = "miRNA", shape = "square",
                 color = list(background = "teal", border="black"))
net <- visGroups(net, groupname = "gene", shape = "dot",       
                 color = list(background = "gold", border="black"))
net
#visLegend(net, main="Legend", position="right", ncol=1) 

# Save plot as PNG
html_name <- tempfile(fileext = ".html")
visSave(net, html_name)

library(webshot); #webshot::install_phantomjs() #in case phantomjs was not installed 

webshot(html_name, zoom = 5, file = "/home/esaha/BONOBO/plots/mRNA_miRNA_breastCancer/LumABsigDiff_network.png")

################################################
# Color genes by KEGG pathway
# Pathways associated to differential genes
pathlist = lapply(pathways[which(names(pathways) %in% sig_pathsAB)], function(x){intersect(x,sig_genes)})
pathlist[which(names(pathlist) %in% sig_path_limma)]
pathlength = lapply(pathlist, length)

# Network Plot
nodes <- data.frame(id = unique(as.vector(as.matrix(diffNet[,c(2,3)]))) , 
                    label=unique(as.vector(as.matrix(diffNet[,c(2,3)]))))
nodes$group <- ifelse(nodes$id %in% diffNet$from, "miRNA", "gene")
node_genes = nodes$id[which(nodes$group == "gene")]

pathtoplot = node_genes
for (i in 1:length(node_genes)){
  g = node_genes[i]
  pathtoplot[i] = names(which.max(pathlength[unlist(lapply(pathlist, function(x){g %in% x}))]))
}
pathtoplot = stringr::str_to_title(lapply(strsplit(pathtoplot, split = "_"), function(x){paste(x[-1], collapse = " ")}))
nodes$group[which(nodes$group == "gene")] = pathtoplot

# Color male edges blue and female edges red
edge_col = rep("red", nrow(diffNet))
edge_col[which(sign(diffNet[,1]) == -1)] = "blue"

diffNet$color = edge_col

net <- visNetwork(nodes, diffNet, width = "100%") 
net <- visGroups(net, groupname = "miRNA", shape = "square",
                 color = list(background = "teal", border="black"))
net
#visLegend(net, main="Legend", position="right", ncol=1) 

# Save plot as PNG
html_name <- tempfile(fileext = ".html")
visSave(net, html_name)

library(webshot); #webshot::install_phantomjs() #in case phantomjs was not installed 

webshot(html_name, zoom = 5, file = "/home/esaha/BONOBO/plots/mRNA_miRNA_breastCancer/LumABsigDiff_network_pathcolor.png")

cbind(node_genes, pathtoplot)

